# }
#
# Update Yold to include all past actions and possibile signings
if(min(Yold$year) > 1950){
# if(is.null(datadir)){
#   datadir <- "~/Dropbox/BiTEN"
# }
# data <- read_env_treaties(datadir, lag=1, write=F, readfile=T)
# Ytemp <- data$Y    # old data
# rm(data)
Ytemp <- load_ytemp(file.path(datadir, "data.RData"))
veryoldrats <- Ytemp[Ytemp$year < min(Yold$year),]
veryoldrats$t <-  veryoldrats$t - max(veryoldrats$t)   # recode t
veryoldrats$j <- Y$j[match(veryoldrats$treaty, Y$treaty)]  # recode j
veryoldrats$i <- Y$i[match(veryoldrats$cowcode, Y$cowcode)]  # recode i
Yold <- rbind(Yold, as.matrix(veryoldrats))
newrats <- Ytemp[Ytemp$year > max(Yold$year) & Ytemp$year <= max(Ynew$year),]   # new ratifications that are not available in current dataset
newrats$t <-  Y$t[match(newrats$year, Y$year)]  # recode t
newrats$j <- Y$j[match(newrats$treaty, Y$treaty)]  # recode j
newrats$i <- Y$i[match(newrats$cowcode, Y$cowcode)]  # recode i
newrats_string <- apply( newrats[, c("cowcode", "j", "t" )], 1, paste0, collapse=",")
oldrats_string <- apply( Ynew[, c("cowcode", "j", "t")], 1, paste0, collapse=",")
suppressWarnings( addrats <- t(sapply(setdiff(newrats_string, oldrats_string), function(z) as.numeric(strsplit(z, ",")[[1]]))) )
keep_add <- apply(addrats[,2:3], 1, function(z) !any(is.na(z)))  # remove j and t indices that are NA
keep_add_string <- apply(addrats[keep_add,], 1, paste0, collapse=",")
newrats1 <- newrats[match(keep_add_string, newrats_string ), ]   # unaccounted for signatures for X
Yold <- rbind(Yold, as.matrix(newrats1))
remove_rows <- unique( which(is.na(Yold[,c("j", "t")]), arr.ind=T)[,1] )   # remove NA rows FOR J AND T
Yold <- Yold[-remove_rows,]
}
# Find initial possible signatures in t0 and those i,j pairs that already signed
possibles <- as.matrix(unique(Ynew[Ynew$t == t0, c("i", "j")]))   # all possible signatures in current year
already_signed <- as.matrix(unique(Yold[Yold$t < t0 & Yold$ratification_year == 1, c("i", "j")]))   # already signed country/treaty pairs
# Find i,j pairs that enter after t0 and the corresponding time when they do so
match_matrix <- matrix(rnorm(S*L), S, L)
match_matrix[rbind( possibles, already_signed)] <- 0   # unsigned and signed country/treaty pairs
leftovers <- which(match_matrix != 0, arr.ind=T)   # not signed AND not possible in year t0
suppressWarnings( tmin <- sapply(1:nrow(leftovers), function(z) min(Ynew$t[Ynew$i == leftovers[z,1] & Ynew$j == leftovers[z,2]]) ) )
tmin[which(is.infinite(tmin))] <- NA   # set infinte to NAs
late_entry <- (tmin > t0 )*tmin
late_entry[which(late_entry <= t0)] <- NA  #NA out non-late entries
late_entry1 <- cbind(leftovers, late_entry)[!is.na(late_entry),]
late_entry <- late_entry1
# Find i,j pairs that exit the dataset early
suppressWarnings( maxt <- sapply(1:nrow(possibles), function(z) max(Ynew$t[Ynew$i == possibles[z,1] & Ynew$j == possibles[z,1]])) )
suppressWarnings( ratt <- sapply(1:nrow(possibles), function(z) Ynew$t[Ynew$ratification_year == 1 & Ynew$i == possibles[z,1] & Ynew$j == possibles[z,1]]) )
maxt[is.infinite(maxt)] <- NA
# ratt[is.infinite(ratt)] <- NA
ratt1 <- rep(NA, length(ratt))
ratt1[which(sapply(ratt, length) > 0)] <- unlist(ratt[which(sapply(ratt, length) > 0)])
early_exit <- cbind(possibles, maxt)[which(maxt < tfinal & is.na(ratt1)), ]  # i,j pairs that exit early
# Both early and late
suppressWarnings( maxt1 <- sapply(1:nrow(late_entry), function(z) max(Ynew$t[Ynew$i == late_entry[z,1] & Ynew$j == late_entry[z,1]])) )
suppressWarnings( ratt1 <- sapply(1:nrow(late_entry), function(z) Ynew$t[Ynew$ratification_year == 1 & Ynew$i == late_entry[z,1] & Ynew$j == late_entry[z,1]]) )
maxt1[is.infinite(maxt1)] <- NA
# ratt1[is.infinite(ratt1)] <- NA
ratt2 <- rep(NA, length(ratt1))
ratt2[which(sapply(ratt1, length) > 0)] <- unlist(ratt1[which(sapply(ratt1, length) > 0)])
both_late_and_early <- cbind(late_entry, maxt1)[which(maxt1 < tfinal & is.na(ratt2)), ]  # i,j pairs that exit early
if(nrow(both_late_and_early) > 0){
both_late_and_early <- cbind( both_late_and_early[,c(1:2)], 0, both_late_and_early[,c(3:4)])
entry_time <- sapply(1:nrow(both_late_and_early), function(z) late_entry[late_entry[,1] == both_late_and_early[z,1] & late_entry[,2] == both_late_and_early[z,2],3])
both_late_and_early[,3] <- entry_time
keep <- apply(both_late_and_early, 1, function(z) !(sum(is.na(z)) > 0))   # remove any NAs
both_late_and_early <- both_late_and_early[keep, ]
early_exit <- rbind(early_exit, both_late_and_early[, c(1,2,4)])  # add to early exit
}
# build possible signatures for all future times and simulations
for(t in trange){
k <- t - t0 + 1
iposs <- cbind(possibles[rep(1:nrow(possibles), times=nsims), ], k, rep(1:nsims, each=nrow(possibles)))   # indices in first year that can be signed
Yhat[,,k,] <- phat[,,k,] <- NA   # all NAs in time period
Yhat[iposs] <- phat[iposs] <- 0   # possible ratifications
jlate_entry <- which(t >= late_entry1[,3])   # AFTER late entry
if(length(jlate_entry) > 0){
# cat("late entry k=", k, "t=", t, "\n")
ilate <- cbind(late_entry1[rep(jlate_entry, times=nsims), 1:2], k, rep(1:nsims, each=length(jlate_entry)) )
Yhat[ilate] <- phat[ilate] <- 0   # 0s for i,j pairs that enter late
}
jearly_exit <- which(t > early_exit[,3])   # AFTER exit
if(length(jearly_exit) > 0){
# cat("early exit k=", k, "t=", t,"\n")
iearly <- cbind(early_exit[rep(jearly_exit, times=nsims), 1:2], k, rep(1:nsims, each=length(jearly_exit)) )
Yhat[iearly] <- phat[iearly] <- NA   # NAs for i,j pairs that leave early
}
}
# # Save all ratifications for future updates
# NAinit1 <- already_signed
# #unique(Yold[Yold$ratification_year == 1, c("i", "j")])   # save already signed treaty/country pairs
# if(!is.null(NApairs)){
#   NAinit1 <- rbind(NAinit1, NApairs)
# }
#
return(list(already_signed=already_signed, possibles=possibles, Yhat=Yhat, phat=phat))
}
r <- pre_predict_build(t0, n, A, B, beta, Yin, Din, X, lag, region_income, seed0=seed0, verbose=FALSE, NApairs=NApairs, filename=outfile, datadir=datadir, outdir=writedir)  # countries=NULL, treaties=NULL,
# Reproduce results from Campbell et. al. (2019) "Latent influence networks in global environmental politics."
#
# Users must:
#   1. Specify parent directory
#   2. Make sure all reuqired libraries in code/libraries.R are installed
#
# Warning: some calls may require large amounts of memory and CPU time to complete
#
rm(list=ls())   # clear current environment
####
# Specify parent directory
# maindir <- "/projects/fmarrs3@colostate.edu/reproduce_envtreat" #
maindir <- "~/Dropbox/BiTEN/reproduce_envtreat"
# "/projects/fmarrs3@colostate.edu/reproduce_envtreat"   # "Specify directory in this string, ~/User/desktop, e.g."
####
#### Source functions
source(file.path(maindir, "code/libraries.R"))
source(file.path(maindir, "code/MLE_functions.R"))
source(file.path(maindir, "code/fit_blin.R"))
source(file.path(maindir, "code/post_functions.R"))
source(file.path(maindir, "code/plot_functions.R"))
####
n <- nsims <- 100
seed0 <- 1
resultsdir="results/final_run"
writedir="results"
cswap=2
tswap=40793
year0=1992
lag = 3
n = nsims
resultsdir <- file.path(maindir, resultsdir)
outdir <- writedir <- file.path(maindir, writedir)
datadir <- file.path(maindir, "data")
if(!dir.exists(writedir)){dir.create(writedir)}
region_income <- read.table(file.path(datadir, "region_income.txt"), header=TRUE)
#### Load data
setwd(resultsdir)
A <- read.table("A.txt")
B <- read.table("B.txt")
beta <- read.table("beta.txt")
load(file.path(datadir, "data.RData"))   # for Y,D,X
load(file.path(resultsdir, "Yhat.RData"))   # for Yhat
####
swap = TRUE
if(swap){
ps <- "_swapped"
t0 <- unique(Y$t[Y$year==year0])+1
NApairs = matrix(c(unique(Y$i[Y$cowcode == cswap]), unique(Y$j[Y$treaty == tswap])), nrow=1)
Din <- D  ;  Yin <- Y
orig_rat <- NULL
for(k in 1:nrow(NApairs)){
orig_rat <- c(orig_rat, Yin$ratification_year[Yin$i == NApairs[k,1] & Yin$j == NApairs[k,2] & Yin$t == (t0-1)])
Din$ratification_year[Din$i == NApairs[k,1] & Din$j == NApairs[k,2] & Din$t == (t0)] <- 0
Yin$ratification_year[Yin$i == NApairs[k,1] & Yin$j == NApairs[k,2] & Yin$t == (t0-1)] <- 0
}
} else {
ps <- "_unswapped"
t0 <- unique(Y$t[Y$year==year0])+1
Din <- D  ;  Yin <- Y
orig_rat <- NULL
NApairs <- NULL
}
# Build old ratification data
pre_predict_build <- function(t0, nsims, A, B, beta, Y, D, X, lag, region_income, tfinal=NULL, response="ratification_year", seed0=1, verbose=F, NApairs=NULL, write_interval=floor(nsims/10), outdir=getwd(), filename=NULL, datadir=NULL)  # countries=NULL, treaties=NULL,
{
if(is.null(tfinal)){ tfinal <- max(Y$t)}
trange <- t0:tfinal
year_range <- unique(Y$year[Y$t == t0] ) : unique(Y$year[Y$t == tfinal] )
# Initialize arrays
Ynew <- Y[Y$t >= t0 & Y$t <= tfinal,]
Yold <- Y[Y$t < t0,]
Dnew <- D[D$t >= t0 & D$t <= tfinal,]
Xnew <- X[X$t >= t0 & X$t <= tfinal,]
S <- max(Dnew$i)  ;   L <- max(Dnew$j)
Yhat <- phat <- array(0, c(S,L,length(trange),nsims))
countries <- sapply(1:S, function(z) unique(Dnew$cowcode[Dnew$i == z]))   # all countries of interest in future
treaties <- sapply(1:L, function(z) unique(Dnew$treaty[Dnew$j == z]))    # all treaties of interest in future
dimnames(Yhat)[[1]] <- dimnames(phat)[[1]] <- countries
dimnames(Yhat)[[2]] <- dimnames(phat)[[2]] <- treaties
dimnames(Yhat)[[3]] <- dimnames(phat)[[3]] <- year_range
region_income$i <- Y$i[match(region_income$cowcode, Y$cowcode)]   # add column for i in region_income indicators
region_income <- region_income[order(region_income$i),]    # reorder
# # Read in big X if X is not big
# if(nrow(X) <  S*L*length(trange)){
#   X <- build_big_X(t0, X, S=S, L=L, tfinal=max(X$t), readfile=T)
# }
#
# Update Yold to include all past actions and possibile signings
if(min(Yold$year) > 1950){
if(is.null(datadir)){
datadir <- "~/Dropbox/BiTEN"
}
data <- read_env_treaties(datadir, lag=1, write=F, readfile=T)
Ytemp <- data$Y    # old data
rm(data)
# Ytemp <- load_ytemp(file.path(datadir, "data.RData"))
veryoldrats <- Ytemp[Ytemp$year < min(Yold$year),]
veryoldrats$t <-  veryoldrats$t - max(veryoldrats$t)   # recode t
veryoldrats$j <- Y$j[match(veryoldrats$treaty, Y$treaty)]  # recode j
veryoldrats$i <- Y$i[match(veryoldrats$cowcode, Y$cowcode)]  # recode i
Yold <- rbind(Yold, as.matrix(veryoldrats))
newrats <- Ytemp[Ytemp$year > max(Yold$year) & Ytemp$year <= max(Ynew$year),]   # new ratifications that are not available in current dataset
newrats$t <-  Y$t[match(newrats$year, Y$year)]  # recode t
newrats$j <- Y$j[match(newrats$treaty, Y$treaty)]  # recode j
newrats$i <- Y$i[match(newrats$cowcode, Y$cowcode)]  # recode i
newrats_string <- apply( newrats[, c("cowcode", "j", "t" )], 1, paste0, collapse=",")
oldrats_string <- apply( Ynew[, c("cowcode", "j", "t")], 1, paste0, collapse=",")
suppressWarnings( addrats <- t(sapply(setdiff(newrats_string, oldrats_string), function(z) as.numeric(strsplit(z, ",")[[1]]))) )
keep_add <- apply(addrats[,2:3], 1, function(z) !any(is.na(z)))  # remove j and t indices that are NA
keep_add_string <- apply(addrats[keep_add,], 1, paste0, collapse=",")
newrats1 <- newrats[match(keep_add_string, newrats_string ), ]   # unaccounted for signatures for X
Yold <- rbind(Yold, as.matrix(newrats1))
remove_rows <- unique( which(is.na(Yold[,c("j", "t")]), arr.ind=T)[,1] )   # remove NA rows FOR J AND T
Yold <- Yold[-remove_rows,]
}
# Find initial possible signatures in t0 and those i,j pairs that already signed
possibles <- as.matrix(unique(Ynew[Ynew$t == t0, c("i", "j")]))   # all possible signatures in current year
already_signed <- as.matrix(unique(Yold[Yold$t < t0 & Yold$ratification_year == 1, c("i", "j")]))   # already signed country/treaty pairs
# Find i,j pairs that enter after t0 and the corresponding time when they do so
match_matrix <- matrix(rnorm(S*L), S, L)
match_matrix[rbind( possibles, already_signed)] <- 0   # unsigned and signed country/treaty pairs
leftovers <- which(match_matrix != 0, arr.ind=T)   # not signed AND not possible in year t0
suppressWarnings( tmin <- sapply(1:nrow(leftovers), function(z) min(Ynew$t[Ynew$i == leftovers[z,1] & Ynew$j == leftovers[z,2]]) ) )
tmin[which(is.infinite(tmin))] <- NA   # set infinte to NAs
late_entry <- (tmin > t0 )*tmin
late_entry[which(late_entry <= t0)] <- NA  #NA out non-late entries
late_entry1 <- cbind(leftovers, late_entry)[!is.na(late_entry),]
late_entry <- late_entry1
# Find i,j pairs that exit the dataset early
suppressWarnings( maxt <- sapply(1:nrow(possibles), function(z) max(Ynew$t[Ynew$i == possibles[z,1] & Ynew$j == possibles[z,1]])) )
suppressWarnings( ratt <- sapply(1:nrow(possibles), function(z) Ynew$t[Ynew$ratification_year == 1 & Ynew$i == possibles[z,1] & Ynew$j == possibles[z,1]]) )
maxt[is.infinite(maxt)] <- NA
# ratt[is.infinite(ratt)] <- NA
ratt1 <- rep(NA, length(ratt))
ratt1[which(sapply(ratt, length) > 0)] <- unlist(ratt[which(sapply(ratt, length) > 0)])
early_exit <- cbind(possibles, maxt)[which(maxt < tfinal & is.na(ratt1)), ]  # i,j pairs that exit early
# Both early and late
suppressWarnings( maxt1 <- sapply(1:nrow(late_entry), function(z) max(Ynew$t[Ynew$i == late_entry[z,1] & Ynew$j == late_entry[z,1]])) )
suppressWarnings( ratt1 <- sapply(1:nrow(late_entry), function(z) Ynew$t[Ynew$ratification_year == 1 & Ynew$i == late_entry[z,1] & Ynew$j == late_entry[z,1]]) )
maxt1[is.infinite(maxt1)] <- NA
# ratt1[is.infinite(ratt1)] <- NA
ratt2 <- rep(NA, length(ratt1))
ratt2[which(sapply(ratt1, length) > 0)] <- unlist(ratt1[which(sapply(ratt1, length) > 0)])
both_late_and_early <- cbind(late_entry, maxt1)[which(maxt1 < tfinal & is.na(ratt2)), ]  # i,j pairs that exit early
if(nrow(both_late_and_early) > 0){
both_late_and_early <- cbind( both_late_and_early[,c(1:2)], 0, both_late_and_early[,c(3:4)])
entry_time <- sapply(1:nrow(both_late_and_early), function(z) late_entry[late_entry[,1] == both_late_and_early[z,1] & late_entry[,2] == both_late_and_early[z,2],3])
both_late_and_early[,3] <- entry_time
keep <- apply(both_late_and_early, 1, function(z) !(sum(is.na(z)) > 0))   # remove any NAs
both_late_and_early <- both_late_and_early[keep, ]
early_exit <- rbind(early_exit, both_late_and_early[, c(1,2,4)])  # add to early exit
}
# build possible signatures for all future times and simulations
for(t in trange){
k <- t - t0 + 1
iposs <- cbind(possibles[rep(1:nrow(possibles), times=nsims), ], k, rep(1:nsims, each=nrow(possibles)))   # indices in first year that can be signed
Yhat[,,k,] <- phat[,,k,] <- NA   # all NAs in time period
Yhat[iposs] <- phat[iposs] <- 0   # possible ratifications
jlate_entry <- which(t >= late_entry1[,3])   # AFTER late entry
if(length(jlate_entry) > 0){
# cat("late entry k=", k, "t=", t, "\n")
ilate <- cbind(late_entry1[rep(jlate_entry, times=nsims), 1:2], k, rep(1:nsims, each=length(jlate_entry)) )
Yhat[ilate] <- phat[ilate] <- 0   # 0s for i,j pairs that enter late
}
jearly_exit <- which(t > early_exit[,3])   # AFTER exit
if(length(jearly_exit) > 0){
# cat("early exit k=", k, "t=", t,"\n")
iearly <- cbind(early_exit[rep(jearly_exit, times=nsims), 1:2], k, rep(1:nsims, each=length(jearly_exit)) )
Yhat[iearly] <- phat[iearly] <- NA   # NAs for i,j pairs that leave early
}
}
# # Save all ratifications for future updates
# NAinit1 <- already_signed
# #unique(Yold[Yold$ratification_year == 1, c("i", "j")])   # save already signed treaty/country pairs
# if(!is.null(NApairs)){
#   NAinit1 <- rbind(NAinit1, NApairs)
# }
#
return(list(already_signed=already_signed, possibles=possibles, Yhat=Yhat, phat=phat))
}
r <- pre_predict_build(t0, n, A, B, beta, Yin, Din, X, lag, region_income, seed0=seed0, verbose=FALSE, NApairs=NApairs, filename=outfile, datadir=datadir, outdir=writedir)  # countries=NULL, treaties=NULL,
#### Source functions
source("~/Dropbox/BiTEN/gitBiTEN/MLE/MLE_functions.R")
# Reproduce results from Campbell et. al. (2019) "Latent influence networks in global environmental politics."
#
# Users must:
#   1. Specify parent directory
#   2. Make sure all reuqired libraries in code/libraries.R are installed
#
# Warning: some calls may require large amounts of memory and CPU time to complete
#
rm(list=ls())   # clear current environment
####
# Specify parent directory
# maindir <- "/projects/fmarrs3@colostate.edu/reproduce_envtreat" #
maindir <- "~/Dropbox/BiTEN/reproduce_envtreat"
# "/projects/fmarrs3@colostate.edu/reproduce_envtreat"   # "Specify directory in this string, ~/User/desktop, e.g."
####
#### Source functions
source("~/Dropbox/BiTEN/gitBiTEN/MLE/MLE_functions.R")
source(file.path(maindir, "code/libraries.R"))
source(file.path(maindir, "code/MLE_functions.R"))
source(file.path(maindir, "code/fit_blin.R"))
source(file.path(maindir, "code/post_functions.R"))
source(file.path(maindir, "code/plot_functions.R"))
####
n <- nsims <- 100
seed0 <- 1
resultsdir="results/final_run"
writedir="results"
cswap=2
tswap=40793
year0=1992
lag = 3
n = nsims
resultsdir <- file.path(maindir, resultsdir)
outdir <- writedir <- file.path(maindir, writedir)
datadir <- file.path(maindir, "data")
if(!dir.exists(writedir)){dir.create(writedir)}
region_income <- read.table(file.path(datadir, "region_income.txt"), header=TRUE)
#### Load data
setwd(resultsdir)
A <- read.table("A.txt")
B <- read.table("B.txt")
beta <- read.table("beta.txt")
load(file.path(datadir, "data.RData"))   # for Y,D,X
load(file.path(resultsdir, "Yhat.RData"))   # for Yhat
####
swap = TRUE
if(swap){
ps <- "_swapped"
t0 <- unique(Y$t[Y$year==year0])+1
NApairs = matrix(c(unique(Y$i[Y$cowcode == cswap]), unique(Y$j[Y$treaty == tswap])), nrow=1)
Din <- D  ;  Yin <- Y
orig_rat <- NULL
for(k in 1:nrow(NApairs)){
orig_rat <- c(orig_rat, Yin$ratification_year[Yin$i == NApairs[k,1] & Yin$j == NApairs[k,2] & Yin$t == (t0-1)])
Din$ratification_year[Din$i == NApairs[k,1] & Din$j == NApairs[k,2] & Din$t == (t0)] <- 0
Yin$ratification_year[Yin$i == NApairs[k,1] & Yin$j == NApairs[k,2] & Yin$t == (t0-1)] <- 0
}
} else {
ps <- "_unswapped"
t0 <- unique(Y$t[Y$year==year0])+1
Din <- D  ;  Yin <- Y
orig_rat <- NULL
NApairs <- NULL
}
# Build old ratification data
pre_predict_build <- function(t0, nsims, A, B, beta, Y, D, X, lag, region_income, tfinal=NULL, response="ratification_year", seed0=1, verbose=F, NApairs=NULL, write_interval=floor(nsims/10), outdir=getwd(), filename=NULL, datadir=NULL)  # countries=NULL, treaties=NULL,
{
if(is.null(tfinal)){ tfinal <- max(Y$t)}
trange <- t0:tfinal
year_range <- unique(Y$year[Y$t == t0] ) : unique(Y$year[Y$t == tfinal] )
# Initialize arrays
Ynew <- Y[Y$t >= t0 & Y$t <= tfinal,]
Yold <- Y[Y$t < t0,]
Dnew <- D[D$t >= t0 & D$t <= tfinal,]
Xnew <- X[X$t >= t0 & X$t <= tfinal,]
S <- max(Dnew$i)  ;   L <- max(Dnew$j)
Yhat <- phat <- array(0, c(S,L,length(trange),nsims))
countries <- sapply(1:S, function(z) unique(Dnew$cowcode[Dnew$i == z]))   # all countries of interest in future
treaties <- sapply(1:L, function(z) unique(Dnew$treaty[Dnew$j == z]))    # all treaties of interest in future
dimnames(Yhat)[[1]] <- dimnames(phat)[[1]] <- countries
dimnames(Yhat)[[2]] <- dimnames(phat)[[2]] <- treaties
dimnames(Yhat)[[3]] <- dimnames(phat)[[3]] <- year_range
region_income$i <- Y$i[match(region_income$cowcode, Y$cowcode)]   # add column for i in region_income indicators
region_income <- region_income[order(region_income$i),]    # reorder
# # Read in big X if X is not big
# if(nrow(X) <  S*L*length(trange)){
#   X <- build_big_X(t0, X, S=S, L=L, tfinal=max(X$t), readfile=T)
# }
#
# Update Yold to include all past actions and possibile signings
if(min(Yold$year) > 1950){
if(is.null(datadir)){
datadir <- "~/Dropbox/BiTEN"
}
data <- read_env_treaties(datadir, lag=1, write=F, readfile=T)
Ytemp <- data$Y    # old data
rm(data)
# Ytemp <- load_ytemp(file.path(datadir, "data.RData"))
veryoldrats <- Ytemp[Ytemp$year < min(Yold$year),]
veryoldrats$t <-  veryoldrats$t - max(veryoldrats$t)   # recode t
veryoldrats$j <- Y$j[match(veryoldrats$treaty, Y$treaty)]  # recode j
veryoldrats$i <- Y$i[match(veryoldrats$cowcode, Y$cowcode)]  # recode i
Yold <- rbind(Yold, as.matrix(veryoldrats))
newrats <- Ytemp[Ytemp$year > max(Yold$year) & Ytemp$year <= max(Ynew$year),]   # new ratifications that are not available in current dataset
newrats$t <-  Y$t[match(newrats$year, Y$year)]  # recode t
newrats$j <- Y$j[match(newrats$treaty, Y$treaty)]  # recode j
newrats$i <- Y$i[match(newrats$cowcode, Y$cowcode)]  # recode i
newrats_string <- apply( newrats[, c("cowcode", "j", "t" )], 1, paste0, collapse=",")
oldrats_string <- apply( Ynew[, c("cowcode", "j", "t")], 1, paste0, collapse=",")
suppressWarnings( addrats <- t(sapply(setdiff(newrats_string, oldrats_string), function(z) as.numeric(strsplit(z, ",")[[1]]))) )
keep_add <- apply(addrats[,2:3], 1, function(z) !any(is.na(z)))  # remove j and t indices that are NA
keep_add_string <- apply(addrats[keep_add,], 1, paste0, collapse=",")
newrats1 <- newrats[match(keep_add_string, newrats_string ), ]   # unaccounted for signatures for X
Yold <- rbind(Yold, as.matrix(newrats1))
remove_rows <- unique( which(is.na(Yold[,c("j", "t")]), arr.ind=T)[,1] )   # remove NA rows FOR J AND T
Yold <- Yold[-remove_rows,]
}
# Find initial possible signatures in t0 and those i,j pairs that already signed
possibles <- as.matrix(unique(Ynew[Ynew$t == t0, c("i", "j")]))   # all possible signatures in current year
already_signed <- as.matrix(unique(Yold[Yold$t < t0 & Yold$ratification_year == 1, c("i", "j")]))   # already signed country/treaty pairs
# Find i,j pairs that enter after t0 and the corresponding time when they do so
match_matrix <- matrix(rnorm(S*L), S, L)
match_matrix[rbind( possibles, already_signed)] <- 0   # unsigned and signed country/treaty pairs
leftovers <- which(match_matrix != 0, arr.ind=T)   # not signed AND not possible in year t0
suppressWarnings( tmin <- sapply(1:nrow(leftovers), function(z) min(Ynew$t[Ynew$i == leftovers[z,1] & Ynew$j == leftovers[z,2]]) ) )
tmin[which(is.infinite(tmin))] <- NA   # set infinte to NAs
late_entry <- (tmin > t0 )*tmin
late_entry[which(late_entry <= t0)] <- NA  #NA out non-late entries
late_entry1 <- cbind(leftovers, late_entry)[!is.na(late_entry),]
late_entry <- late_entry1
# Find i,j pairs that exit the dataset early
suppressWarnings( maxt <- sapply(1:nrow(possibles), function(z) max(Ynew$t[Ynew$i == possibles[z,1] & Ynew$j == possibles[z,1]])) )
suppressWarnings( ratt <- sapply(1:nrow(possibles), function(z) Ynew$t[Ynew$ratification_year == 1 & Ynew$i == possibles[z,1] & Ynew$j == possibles[z,1]]) )
maxt[is.infinite(maxt)] <- NA
# ratt[is.infinite(ratt)] <- NA
ratt1 <- rep(NA, length(ratt))
ratt1[which(sapply(ratt, length) > 0)] <- unlist(ratt[which(sapply(ratt, length) > 0)])
early_exit <- cbind(possibles, maxt)[which(maxt < tfinal & is.na(ratt1)), ]  # i,j pairs that exit early
# Both early and late
suppressWarnings( maxt1 <- sapply(1:nrow(late_entry), function(z) max(Ynew$t[Ynew$i == late_entry[z,1] & Ynew$j == late_entry[z,1]])) )
suppressWarnings( ratt1 <- sapply(1:nrow(late_entry), function(z) Ynew$t[Ynew$ratification_year == 1 & Ynew$i == late_entry[z,1] & Ynew$j == late_entry[z,1]]) )
maxt1[is.infinite(maxt1)] <- NA
# ratt1[is.infinite(ratt1)] <- NA
ratt2 <- rep(NA, length(ratt1))
ratt2[which(sapply(ratt1, length) > 0)] <- unlist(ratt1[which(sapply(ratt1, length) > 0)])
both_late_and_early <- cbind(late_entry, maxt1)[which(maxt1 < tfinal & is.na(ratt2)), ]  # i,j pairs that exit early
if(nrow(both_late_and_early) > 0){
both_late_and_early <- cbind( both_late_and_early[,c(1:2)], 0, both_late_and_early[,c(3:4)])
entry_time <- sapply(1:nrow(both_late_and_early), function(z) late_entry[late_entry[,1] == both_late_and_early[z,1] & late_entry[,2] == both_late_and_early[z,2],3])
both_late_and_early[,3] <- entry_time
keep <- apply(both_late_and_early, 1, function(z) !(sum(is.na(z)) > 0))   # remove any NAs
both_late_and_early <- both_late_and_early[keep, ]
early_exit <- rbind(early_exit, both_late_and_early[, c(1,2,4)])  # add to early exit
}
# build possible signatures for all future times and simulations
for(t in trange){
k <- t - t0 + 1
iposs <- cbind(possibles[rep(1:nrow(possibles), times=nsims), ], k, rep(1:nsims, each=nrow(possibles)))   # indices in first year that can be signed
Yhat[,,k,] <- phat[,,k,] <- NA   # all NAs in time period
Yhat[iposs] <- phat[iposs] <- 0   # possible ratifications
jlate_entry <- which(t >= late_entry1[,3])   # AFTER late entry
if(length(jlate_entry) > 0){
# cat("late entry k=", k, "t=", t, "\n")
ilate <- cbind(late_entry1[rep(jlate_entry, times=nsims), 1:2], k, rep(1:nsims, each=length(jlate_entry)) )
Yhat[ilate] <- phat[ilate] <- 0   # 0s for i,j pairs that enter late
}
jearly_exit <- which(t > early_exit[,3])   # AFTER exit
if(length(jearly_exit) > 0){
# cat("early exit k=", k, "t=", t,"\n")
iearly <- cbind(early_exit[rep(jearly_exit, times=nsims), 1:2], k, rep(1:nsims, each=length(jearly_exit)) )
Yhat[iearly] <- phat[iearly] <- NA   # NAs for i,j pairs that leave early
}
}
# # Save all ratifications for future updates
# NAinit1 <- already_signed
# #unique(Yold[Yold$ratification_year == 1, c("i", "j")])   # save already signed treaty/country pairs
# if(!is.null(NApairs)){
#   NAinit1 <- rbind(NAinit1, NApairs)
# }
#
return(list(already_signed=already_signed, possibles=possibles, Yhat=Yhat, phat=phat))
}
r <- pre_predict_build(t0, n, A, B, beta, Yin, Din, X, lag, region_income, seed0=seed0, verbose=FALSE, NApairs=NApairs, filename=outfile, datadir=datadir, outdir=writedir)  # countries=NULL, treaties=NULL,
r <- pre_predict_build(t0, n, A, B, beta, Yin, Din, X, lag, region_income,
seed0=seed0, verbose=FALSE, NApairs=NApairs,
filename=outfile, datadir=NULL, outdir=writedir)  # countries=NULL, treaties=NULL,
# datadir
already_signed <- r$already_signed
possibles <- r$possibles
Yhat <- r$Yhat
phat <- r$phat
save(already_signed, possibles, Yhat, phat, file=file.path(datadir, "pre-prediction2.RData"))
rm(list=ls())
datadir <- "~/Dropbox/BiTEN/reproduce_envtreat/data"
setwd(datadir)
load("pre-prediction.RData")
already_signed1 <- already_signed
early_exit1 <- early_exit
late_entry11 <- late_entry1
possibles1 <- possibles
load("pre-prediction2.RData")
range(already_signed1 - already_signed)
range(already_signed1 - already_signed, na.rm=TRUE)
range(early_exit1 - early_Exit, na.rm=TRUE)
range(early_exit1 - early_exit, na.rm=TRUE)
range(late_entry11 - late_entry1, na.rm=TRUE)
early_exit
early_exit1
range(possibles1 - possibles, na.rm=TRUE)
# Reproduce results from Campbell et. al. (2019) "Latent influence networks in global environmental politics."
#
# Users must:
#   1. Specify parent directory
#   2. Make sure all reuqired libraries in code/libraries.R are installed
#
# Warning: some calls may require large amounts of memory and CPU time to complete
#
rm(list=ls())   # clear current environment
####
# Specify parent directory
# maindir <- "/projects/fmarrs3@colostate.edu/reproduce_envtreat" #
maindir <- "~/Dropbox/BiTEN/reproduce_envtreat"
# "/projects/fmarrs3@colostate.edu/reproduce_envtreat"   # "Specify directory in this string, ~/User/desktop, e.g."
####
#### Source functions
source("~/Dropbox/BiTEN/gitBiTEN/MLE/MLE_functions.R")
source(file.path(maindir, "code/libraries.R"))
source(file.path(maindir, "code/MLE_functions.R"))
source(file.path(maindir, "code/fit_blin.R"))
source(file.path(maindir, "code/post_functions.R"))
source(file.path(maindir, "code/plot_functions.R"))
####
